このページの目的

慶應SFC 2017年秋学期開口科目Data Science for Environmental Governance(Spatial Data Modeling)においてあまり扱わなかった空間データとcsvデータを結合し、空間データ分析を行う方法を解説することです。 具体的にはRを用いてshapeファイルとcsvデータを結合し、図に表し、Moran’s Iを計算するぐらいまでは解説したいと思います。

利用するデータ

今回利用するデータはGitHubにて公開しています(https://github.com/N-Yukihiro/lectute_2017_spatial)。 shapeファイルは国土数値情報ダウンロードサービス(http://nlftp.mlit.go.jp/ksj/index.html)より東京都の行政区域をダウンロードしました。 そして、今回は空間データ分析によく利用される地価を分析したいと思います。 地価は平成29年地価公示価格(東京都分)(http://www.zaimu.metro.tokyo.jp/kijunchi/29kouji/index.html)より23区に該当する地域の住宅地価をcsvに直してlandprice.csvとして保存しています。

実践

パッケージのインストールと読み込み

まずは今回利用するパッケージをインストールし、読み込みます。

install.packages(c("readr", "dplyr", "sf", "devtools","methods", "spdep"))
devtools::install_github("tidyverse/ggplot2")
library(readr)
library(dplyr)
library(sf)
library(ggplot2)
library(methods)
library(spdep)

readrはcsvなどを読み込むため、dplyrはデータフレームを操作するため、sfは空間データを操作するため、ggplot2は空間データを図に表すため、methodssfパッケージで読み込んだデータをspdepパッケージで扱えるように変換するため、spdepは空間データ分析のために使います。 なお、ggplot2は開発版を使います。

データの読み込み

land_price <- readr::read_csv("landprice.csv", locale = locale(encoding = "cp932")) %>% 
  mutate(code_f = as.factor(code))

上のコードではreadrパッケージの中のread_csv関数を利用してcsvファイルを読み込みland_priceという変数に格納しています。 locale=locale(encoding="cp932")でShift-JISで読み込むことを指定しています。 次の行で整数値として格納されている市町村コードをファクターに直しています。

データとしてはこのようになっています。

Dist <- sf::st_read(dsn = "shape", layer = "N03-17_13_170101",options = "ENCODING=cp932") %>% 
  dplyr::group_by(N03_007) %>%
  dplyr::summarise(geometry = st_union(geometry)) %>%
  dplyr::ungroup()

このコードではsfパッケージの中のst_read関数を利用してshapeファイルを読み込みDistという変数に格納しています。 dsnでshapeファイルが入っているフォルダを指定し、layerでshaperファイルを指定します。 Rのコードとshapeファイルが同じフォルダに入っている場合にはdsn = getwd()とします。 layerを指定するときに拡張子は書きません。

2行名以降では海を跨いでマルチフィーチャーになっているデータをシングルフィーチャーに直しています。

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## http://rstudio.github.io/DT/server.html

データの部分(dbfファイルにあたる)はこのようになっています。

データの結合

match <- dplyr::inner_join(Dist, land_price, by = c("N03_007" = "code_f"))

読み込んだcsvとshapeファイルを1つのファイルに結合する操作をしています。 そのためにdplyrパッケージの中のinner_join関数を利用しています。 はじめの2つの引数でマッチさせたいデータを指定します。 最初の引数にshapeファイルを指定して、第二引数にcsvファイルを指定しましょう。

そして、inner_join関数を使う際にキーになる列名を指定します。 今回は列名がそれぞれのデータで異なるので、by = c("column_name1" = "column_name2")とします。

inner_join関数では、お互いのデータセットで共通するキーが存在するもののみデータを残します。 詳しくはdplyrの使い方を調べてみてください。

データの表示

match %>% 
  dplyr::select(LandPrice) %>% 
  plot()

まずはオーソドックスなplotによる表示です。 plotでは全ての列を表示してしまうので、selectを使って表示する列だけを選択しています。

plot(match["LandPrice"])

もっと簡単に上のようなコードで表示することも可能です。

ggplot2::ggplot(data = match) +
  ggplot2::geom_sf(aes(fill = LandPrice))

最後にggplot2による地図の表示の仕方を紹介します。 sfパッケージで読み込んだファイルはggplot2パッケージのgeom_sf関数で簡単に表示することができます。 まず、ggplot(data=data)で利用するデータを指定し、geom_sf(aes(fill=column_name))で表示したい列名を指定します。 geom_sfでは、spパッケージなどのsfパッケージ以外で読み込んだファイルは表示できないことに注意してください。

私はtmapによる描画もけっこう好きなのですが、なぜか私のWindowsで動かないので割愛します。 ちなみに普段はBash on Ubuntu on Windowsでtmapを動かしています。

Moran’s Iの計算

空間的自己相関を表す代表的な指標であるMoran’s Iを計算してみましょう。 空間的自己相関を計算するためにはまず接続行列を与える必要があります。

coords <- sp::coordinates(methods::as(match, "Spatial"))
match_Delaunay_nb <- spdep::tri2nb(coords)
match_simple_nb <- spdep::poly2nb(methods::as(match, "Spatial"), queen=F)

まず、sfパッケージで読み込んだデータをspdepが扱えるようにデータを変換してmatchファイルの緯度経度を求めます。 そして、tri2nbを使ってドロネーの三角網によって接続関係を与えます。 poly2nbを使うと単純接続によって隣接を与えます。

plot(methods::as(match, "Spatial"))
par(new=T)
plot(match_Delaunay_nb,coords=coords, add=TRUE, pch=16, col='darkred')

ドロネーの三角網によって、このように接続関係が与えられています。

plot(methods::as(match, "Spatial"))
par(new=T)
plot(match_simple_nb,coords=coords, add=TRUE, pch=16, col='darkred')

単純接続によって、このように接続関係が与えられています。

それではこれらを使ってMoran’s Iを求めてみましょう。 Moran’s Iは次の数式の通りです。

\[ I = \frac{n}{S_0}\frac{\Sigma_{i=1}^{n}\Sigma_{j=1}^{n}w_{ij}(y_i-\overline{y})(y_j-\overline{y})}{\Sigma_{i=1}^{n}(y_i-\overline{y})}. \]

spdep::moran.test(match$LandPrice, nb2listw (match_Delaunay_nb, style="W"))
## 
##  Moran I test under randomisation
## 
## data:  match$LandPrice  
## weights: nb2listw(match_Delaunay_nb, style = "W")  
## 
## Moran I statistic standard deviate = 4.285, p-value = 9.135e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.377998162      -0.045454545       0.009765631

ドロネーの三角網によって与えた接続を利用してMoran’s Iを計算します。

spdep::moran.test(match$LandPrice, nb2listw (match_simple_nb, style="W"))
## 
##  Moran I test under randomisation
## 
## data:  match$LandPrice  
## weights: nb2listw(match_simple_nb, style = "W")  
## 
## Moran I statistic standard deviate = 3.7406, p-value = 9.179e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.33950161       -0.04545455        0.01059110

最後に単純接続によって与えた接続を利用してMoran’s Iを計算します。

以上で簡単でしたが、shapeファイルとcsvを結合して空間データ分析をするための解説でした。 質問があれば、ご連絡ください。

その他のRの使い方について

Rの使用方法については日本語、英語を問わず多くの資料がインターネット上に公開されています。

一応今までに私が作成したスライドなどはhttps://n-yukihiro.github.io/post/slides/にまとめています。 具体的にはdplyrggplot2の使い方などを示しています。

Enjoy!